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We continue our exploration of whether the flyby anomahes can be explained by scattering 
of spacecraft nucleons from dark matter gravitationally bound to the earth. We formulate 
and analyze a simple model in which inelastic and elastic scatterers populate shells generated 
, by the precession of circular orbits with normals tilted with respect to the earth's axis. Good 

• fits to the data published by Anderson et al. are obtained. 
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I. INTRODUCTION 

CM 

r~| ' In this paper we follow up our earlier investigation [1] of the anomalous geocentric frame orbital 



energy changes that are observed during earth flybys of various spacecraft, as reported by Anderson 
et al. Some flybys show energy decreases, and others energy increases, with the largest 

anomalous velocity changes of order 1 part in 10^. While the possibility that these anomalies are 
. artifacts of the orbital fitting method used in is still being actively explored, there is also a chance 

' that they may represent new physics. In {l| we explored the possibility that the flyby anomalies 



T^lj- . result from scattering of spacecraft nucleons from dark matter particles in orbit around the earth, 

(N ■ 

with the observed velocity decreases arising from elastic scattering, and the observed velocity 



increases arising from exothermic inelastic scattering, which can impart an energy impulse to a 
spacecraft nucleon. Many constraints on this hypothesis were analyzed in [1], with the conclusion 
that the dark matter scenario is not currently ruled out, but requires dark matter to be non-self- 



00 
O 

o\ 
o 

> 

X 

H annihilating, with the dark matter scattering cross section on nucleons much larger, and the dark 

matter mass much lighter, than usually assumed. 

However, no attempt was made in jl| to construct a model for the spatial and velocity distribu- 
tion functions for dark matter populations in earth orbit, to see whether it can fit the flyby data 
reported in 3l- Formulating such a model is the aim of the present paper. Our basic assumption is 
to consider two populations of dark matter particles, one of which scatters on nucleons elastically, 
and the other of which scatters inelastically, each with a shell-like distribution of orbits generated 
by the precession of a tilted circular orbit around the earth's rotation axis. The formulas defining 
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this model are developed in Sec. II, with details of derivations in Appendices, and the results of 
numerical fits to the fiyby data are given in Sec. III. We show that good fits to the data are possi- 
ble, which leaves dark matter scattering as a viable candidate for explaining the flyby anomalies, 
pending further investigation of possible artifactual explanations^ of the flyby data, and further 
experiments aimed at directly detecting dark matter and determining its properties. 

II. FORMULAS DEFINING THE MODEL 
A. Velocity change formulas 

n 

We recall from [1] formulas for the velocity change when a spacecraft nucleon of mass nii ~ 
IGeV and initial velocity ui scatters from a primary dark matter particle of mass m2 and initial 
velocity U2, into an outgoing nucleon of mass mi and velocity vi, and an outgoing secondary 
dark matter particle of mass m'2 = m2 — Am and velocity V2 ■ The inelastic case corresponds 
to 7^ m2, while in the elastic case, = m2 and Am = 0. Under the assumptions, (i) both 
initial particles are nonrelativistic, so that |ni| << c, |n2| << c, (ii) the center of mass scattering 
amplitude f{9) depends only on the auxiliary polar angle 9 of scattering, and (iii) in the exothermic 
inelastic case, Am/m2 and 7712/7712 are both of order unity, a straightforward calculation gives the 
outgoing nucleon velocity change, averaged over scattering angles. In the elastic scattering case, 
with Am = 0, m2 = m2, we have 

{dvi) = -2 (ni-^2)(sin^(g/2)) , (1) 

mi + 7n,2 

while in the inelastic case a good approximation is 

(\ 1/2 
2A 777 m'2 \ 
-. 7- c{cos9} , (2) 
mi (mi + mgj / 

with (...) denoting the angular average over the center of mass differential scattering cross section. 
Since ui and U2 are typically of order 10 km s~^, the velocity change in the inelastic case is 
significantly larger than that in the elastic case. 



^ A parameterized post-Newtonian analysis, given in an unpublished memo in the "Talks+Memos" section of the 
author's home page, shows that deviations from Einstein gravity within the framework of metric theories of gravity 
obeying the equivalence principle cannot give residual accelerations large enough to explain the flyby anomalies. 
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B. Change in outgoing spacecraft velocity 

Again as shown in to get the force per unit spacecraft mass resulting from dark matter 
scatters, that is, the acceleration, one multiplies the velocity change in a single scatter {6vi) by the 
number of scatters per unit time. This latter is given by the flux |ui — U2\, times the scattering 
cross section a, times the dark matter spatial and velocity distribution p(^x,U2)- Integrating out 
the dark matter velocity, one thus gets for the force acting at the point x{t) on the spacecraft 
trajectory with velocity ui = dx{t)/dt, 



SF = J d^U2{6vi)\ui - U2\a-p{x,U2) ■ (3) 

Equating the work per unit spacecraft mass along a trajectory from ti to tj to the change in kinetic 
energy per unit mass (assuming that the initial and final times are in the asymptotic region where 
the potential energy can be neglected) we get 



6-{vf - v1) =vf • 5vf = I dt{dx/dt) ■ 6F 



dt J d U2{dx/dt) ■ {6vi)\ui — U2\o'p{x,U2) 

(4) 

C. Cross section and scattering-angle averaged kinematics 

Let W be the center of mass scattering energy of the dark matter-spacecraft nucleon system. 
A simple calculation shows that to a good approximation we have 

W ^ ^ mim2 {ui - U2Y , . 



(ttt-i + 1712)0^ 2(mi + 7712)^ 

and so for 7712 < mi and for the nonrelativistic velocities ui, U2 of interest, the scattering is very 
close to threshold. Thus the cross section will be dominated by the lowest partial waves, which near 
threshold each have a characteristic power law dependence on the entrance channel momentum 

\Ui - U2\ . (6) 



nil + 

For elastic scattering, the cross section is S-wave dominated, and tends to a fe-independent 
constant a^i near threshold, and the angular average 2(sin2(6l/2)) reduces to 1 - (cos 6*) = 1. Thus 
when Eq. ([1]) is substituted into Eqs. 1^ and ([!]), we can effectively replace 2(sin^(0/2))(7 by the 
/c-independent constant iTei- 
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For exothermic inelastic scattering, the leading contribution to {cos 9) comes from the inter- 
ference term between the S- and P-waves in the cross section, which scales [2] as k^'^k^^'^k^/'^ ~ 
constant near threshold. Writing near threshold 

we have 

(COS0) ~Ancl/(AnclA:"^) ■ 

(8) 

So when Eq. ^ for the inelastic exothermic case is substituted into Eqs. ([3]) and ([4]), we can 
effectively replace {cos 9)a by the /c-independent constant -Bineh remembering, however, that this 
is not the total cross section (which approaches ^inei^~^ near threshold) but is proportional to the 
coefficient of the S-wave P-wave interference term in the differential cross section. 



D. The dark matter distribution function p{x, U2) 

We now address the task of formulating a model for the distribution function p{x,U2) that 
describes dark matter postulated to be in orbit around the earth. The simplest model would be a 
disk composed of dark matter in circular orbits in earth's equatorial plane, but attempts to fit the 
flyby anomaly data with such a model were unsuccessful, since for any reasonable disk inner radius, 
some of the flybys (such as NEAR) pass inside the disk. We thus proceed to the next simplest 
model, which is constructed from dark matter in a circular orbit, of radius r and tilted at an angle 
tp (0 < ip < tt) with respect to earth's equatorial plane. If the earth were exactly spherically 
symmetric, its gravitational field would be strictly monopole, and such a tilted orbit would be 
stable. But in fact the earth's rotation produces an equatorial bulge, and so its mass distribution 
is only axially symmetric around its rotation axis, giving rise to quadrupole and higher moments 
in its gravitational field. As a result of these higher moments, the tilted orbit precesses around the 
earth's rotation axis, in such a way that the angular momentum component along the earths's 
axis is conserved. Over a long period of time, this precession will smear an initial cluster of tilted 
orbits into a uniform shell, obtained by averaging the tilted circle over the azimuthal angle that 
its normal makes with respect to the earth's rotation axis. 

To give this picture a mathematical description, let x, y,z he a Cartesian axis system, with 
positive z pointing to the earth's North pole (so that the rotation sense of the earth is from 
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X to y). Let the normal h to the tilted orbit have polar angle ■0 and azimuthal angle (j) with 
respect to this system, so that n{ip, (j)) = (sin ifj cos sin tjj sin cos '(/'), and let the angle of rotation 
within the plane of the dark matter orbit be 9, with increasing 9 corresponding, at = 0, to the 
direction of earth's rotation. Then a parametric description of the tilted circle is P{r,9,<j3) = 
{P^{r, 9, ^),Py{r, 9, <P),P, (r, 9, <P)) , with 

Px{r, 9, (j)) = r(cos 9 cos V' cos </> — sin sin (j)) , 
Py{r,9,(j)) = r {cos 9 cos ip sin (j) + sin 9 cos (j)) , 
Pz{r, 9, <j)) = — r cos 9 sin ip , \Pz{r, 9, (j))\ < r sin ip 

(9) 

The corresponding velocity unit vector of a dark matter particle in the tilted circular orbit is 

U{9,cj)) = {U,{9,(l)),Uy{9,cl)),U,{9,cl))) =r-^dP/d9, with 

Ux = — sin 9 cos ip cos (j) — cos 9 sin (j) , 
Uy = — sin 9 cos ip sin (p + cos 9 cos (j) , 

U^= sin 61 sin • (10) 

The velocity vector is obtained by multiplying the velocity unit vector by the velocity magnitude 
(GM©/r)^/2 fQj. a particle in a circular orbit of radius r, with G the Newton gravitational constant 
and the earth mass. 

Integrating the position and velocity distribution for a tilted circular orbit over the angles 9, (j) 
gives the distribution for the corresponding shell, and integrating over the shell parameters r, tp 
with a general weighting function w{r, ip) gives as the model for the dark matter distribution 
function 

dr j dipw{r,ip) j d9 j d(P5^{x-P{r,9,<p))5^{u2-{GM(Q/rf/^U{9,(P)) , (11) 
with the corresponding total number of particles in the shell given by 

A^ = j j d^U2p{x,U2) =4^^y" dr j dipw{r,ip) . (12) 
Referring to Eq. dH), we have to evaluate an integral over the distribution function of the form 

/ = j dt j d^U2F{x{t),dx{t)/dt,U2)p{x{t),U2) , (13) 
with F{x{t),dx{t)/dt,U2) given by 

F{x{t),dx{t)/dt,U2) = {dx{t)/dt) ■ {6vi)\u^=dx{t) / dt \dx{t)/dt - U2\a . (14) 
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On substituting Eq. (jllh and noting that the coordinate delta function constrains r = \P{r,6, 
\x{t)\ = r{t), we obtain 



dt J dip w{r{t),ip) J dr 

f2n f-2n 

X del d(pF{x{t),dx{t)/dt,{GMs)/r{t))^/^U{9,(l)))6''{x{t) - P{r,e,^)) 
Jo Jo 



(15) 



As shown in Appendix A, by making changes of variable one can carry out the integrations over 
r, (j) and 9 in Eq. (|15p. leaving an integral in which 9 and z have been replaced, by virtue of the 
delta function constraints, by 6{x{t)) and z{t) = z{x(t)), 

1 = j dt j di;w{r{t),i;)Y,F{m,dx{t)/dt,{GM^/r{t))^/^U±{9{^^^^^ 
1 



sin^ V — 

(16) 

Note that by virtue of Eq. ([9]), the integration domain extends only over \z{t)\ < r{t)smip, and 
hence the argument of the square root is nonnegative. In Eq. (I16p the sum over it is over the two 
roots 9{x(t)) of the equation cos9{x{t)) = — z(t)/(r(t) sin -0), which differ in the sign of sm9, 

sm9{x{t)) = ±^Jl- z{tf/{r{t)^sm'^ij) , (17) 

while the values of (j){x{t)) corresponding to these two roots 9{x{t)) are obtained by equating x{t) 
to P and then solving Eq. ^ for cos(/» and sin(/>. The two roots correspond to the fact that a 
circular orbit with tilt angle ip consists of two semicircular segments, with opposite directions of 
the velocity component normal to the equatorial plane. Thus the intersection of the spacecraft 
trajectory x{t) with the dark matter shell generated by azimuthal rotation of such a tilted circular 
orbit will intersect two segments of circular orbits, one up-going and one down-going relative to 
the equatorial plane. 

It will be useful for what follows to express the unit velocities U±{9{x{t)) , (j){x{t))) in terms of 
their components on unit vectors n^\{t) = (z x x{t))/\z x x{t)\ and 'h±{t) = x{t) x h^\{t), normal to 
x{t) = x{t)/r, that are respectively parallel (in the sense of earth rotation) and perpendicular to 
the earth equatorial plane. A simple calculation given in Appendix B shows that U± are given on 
this basis by 

U±i9{x{t)), ^{x{t))) = C{t)hu ± D{t)hi_ , (18) 
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with the coefficients C{t) and D{t) given by 



^ ^ rjt) cosj' ^ ^ v^r(t)2sin^V>- 

which obey C(t)2 + D{tf' = 1. Explicit expressions for n\\{t) and ri_L(t) in the flyby plane basis 
are given in the next subsection, which together with Eqs. (jlSp and (jl9p give the formulas for the 
unit velocities U±i9{x{t))^ (t){x{t))) on the flyby plane basis needed in the numerical computations. 



E. Flyby orbital plane kinematics 

The Anderson et al. paper Q] gives the flyby orbit parameters in terms of coordinates on 
the celestial sphere, but it will be more convenient for our purposes to carry out all flyby orbit 
calculations in the flyby orbital plane. Let Xq, Vo^ Zq be a Cartesian axis system, with Zq normal to 
the flyby orbital plane. The flyby orbit can then be written in parametric form as 

Xa{t) =r{t) cos 9 o{t) , 
yo{t)=r{t) sin 9o{t) , 



l + ecosOoit) ^ 1 + e 

^ m/^^ -V/sin6lo(t) -yo{t) ^a(.\Ku 

= 1 + e = l + ecosO^itf ^^'^'''' ' 

1 + e l + ecost^o(i) 

RfVf 



(20) 



The scale parameter p, the eccentricity e, the velocity at closest approach to earth V/, the radius 
at closest approach^ii;, and the velocity at inftai.y a.c given in Table I for each of the six 

flybys discussed in [2], together with the polar angle / and azimuthal angle a of the earth's north 
pole with respect to the Xo,yo,Zo coordinate system. The quantities Vf and Voo are given directly 



while Rf,p, and e can be calculated from them using the formulas 



e=l + 



2K3 



4G'Ma 

oo 



* rv 



+ 



* rv 



V? - T/2 



The earth axis polar angle / is also directly given in 
calculated from the formula 



(21) 



while the azimuthal angle a can be 



sin{ 



cos a 



sin J 



(22) 



with (j)' the geocentric latitude at closest approach (which is called in [3]; with the orbit 
parametrization of Eq. ()20p . cp' is the latitude of the positive Xq axis). This formula does not 
determine the c^uadrant in which a lies, but this can be fixed from the additional orbital param- 
eters given in [2:] (with some corrections supplied to me by J.K. Campbell [4]). Enough orbit 
parameters are given in \^ to provide several redundancies that serve as cross-checks on these 
calculations. 



TABLE I: Flyby orbital parameters 





GLL-I 


GLL-II 


NEAR 


Cassini 


Rosetta 


Messenger 


Vf (km/s) 


13.740 


14.080 


12.739 


19.026 


10.517 


10.389 


Rf (km) 


7,334 


6,674 


6,911 


7,544 


8,332 


8,715 


Voo (km/s) 


8.949 


8.877 


6.851 


16.010 


3.863 


4.056 


e 


2.474 


2.320 


1.814 


5.851 


1.312 


1.360 


p (km) 


25,480 


22,160 


19,450 


51,690 


19,260 


20,570 


/ (deg) 


142.9 


138.7 


108.0 


25.4 


144.9 


133.1 


a (deg) 


-45.1 


-147.4 


-55.1 


-158.4 


-53.1 


0.0 
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To carry out the computation of the fly by velocity change in the flyby plane basis Xo,yo, Zq we 
will need the components of ny and n± on this basis. In Eq. ()Bip we gave their components on 
the earth centered basis x, y, z; these can be rotated to the flyby plane basis, but it is simpler to 
calculate them directly by going back to the defining cross product relations, using the components 
of x{t) and of the earth axis z on the flyby plane basis, 



x{t) ={xo{t),yoit),0) 



z =(sin I cos Q, sin / sin q, cos /) 



(23) 



From these we find 




( ~ yo{t) cos /, Xo{t) cos /, {yo{i) cos a — Xo{t) sin q) sin /) 



X {yo{t){yo{t) cos a — Xo{t) sin a) sin I, —x, 



o 



{t){yo{t) cos a — Xo{t) sin a) sin I, r(t)^ cos /) 



(24) 



with 



r{t) =\x{t)\ = V^oitr+Voit? 



z{t) =x(t) ■ z = {xo{t) COS a + yo{t) sin a) sin / 



(25) 



Substituting Eq. ([20]) for Xo{t) and yo{t) into Eq. ([25]) we have 



z{t) = r{t) sin / cos {6o{t) — a 



) 



(26) 



which allows one to rewrite the Jacobian factor appearing in Eq. (|16p as 



1 1 



(27) 




when the argument of the square root is nonnegative. 



F. Simplified model used for numerical work 



The model as defined above involves a general weighting function w{r,ij))^ but for an initial 
survey we make the simplifying assumption of only a single tilt angle ipi, ^ip^ for the inelastic and 



10 



elastic scatterers, respectively, and Gaussian distributions in r with different centers and widths 
for each. Thus we take for the inelastic scatterers 

Wi{r,^) = Kie-^^-''^^"/''h{^l;-^i) , (28) 

and for the elastic scatterers 

We{r, ij) = iv:ee-(^-^=)'/^'(5(V' - V-e) • (29) 

With this choice, the integral of Eq. (jl2p becomes 

Ni = A-K^/^KeDi , £ = i,e . (30) 

It is now convenient to combine the constants Ki^^ with the mass-dependent constants appearing 
in Eqs. ([1]) and ([2|) of Sec. IIA, and the constants a^i and Anei introduced in Sec. IIC, giving 
new parameters pi, pe characterizing the effective density times cross section for the inelastic and 
elastic scatterer distributions, 



Pe = ; CTcl-Ke , 

mi + 1712 



2Am m'r 



1/2 



Pi = 7 ; 7T B^^c\Ki . 

\ mi (mi + mgj / 

(31) 

Thus in Eq. ([15]) we effectively replace (see Eq. ([HI)) 

j d'ipw{r{t),'ilj) F{x{t),dx{t)/dt,U2) (32) 

by 

J2 {\dxit)/dt-U2\{dx{t)/dt) ■ l^£/o^e-('-W-^^)'/^l}|^=^, , (33) 

£=i,e 

with Vi given by 

Vi =c [dx{t)/dt — U2) /\dx{t)/dt — U2\ , 
Ve = - {dx{t) / dt - U2) , 

(34) 

and with U2 evaluated as U± of Eqs. ([16|) and ([T8|) . The simplified model thus defined has 
eight parameters, four parameters V'i) Pii Rii Di characterizing the inelastic scatterers, and four 
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parameters Po Re, characterizing the elastic scatterers. Finahy, we note that by combining 
Eqs. (fSOl) and (fSTI) . and approximating 

1/2 



1712 I 2Am m'2 \ 771-2 



mi + 1712 \ TTii (mi + 771-2) / "^1 ' 

we find the fohowing estimates for the total mass in the dark matter shells, 



(35) 



(36) 



III. NUMERICAL RESULTS AND DISCUSSION 



Let us turn now to numerical fitting of the eight parameter model to the flyby anomalies reported 
in Qj. In carrying out the needed integrals over flyby orbits, we replaced the integration over t by 
an integration over orbit angle 9o, using the expression for dOo/dt given in Eq. (I20p . To utilize 
integration mesh points efficiently, the integrations were restricted to the parts of the orbits where 
the Gaussian factors e~*-^~^*^^/^^ were larger than = 0.00012, that is, to the parts of the orbits 
where |r — Ri\ < 2>D^. 

In attempting to search for good fits with coarse meshes, we found that the infinite jump in the 
Jacobian factor of Eq. (I27p at the dark matter shell edges led to the search program settling on 
false minima reflecting truncation errors, which were unstable with respect to small changes in the 
integration mesh or fitting parameters. To avoid this problem, we replaced the original Jacobian 
by a smoothed Jacobian, as follows. Abbreviating W = z (t)"^ /{r{t)^ sin"^ tp), and using to denote 
the usual step function, the original Jacobian contains the function with an infinite jump atW = l, 

nW) = ?^ . (37) 

We replaced this by the following function, which is continuous and has a continuous first derivative, 

/e(VF) =^e-^^(^) for > 1 - e , 



= - ^(^ - 1 + e) + ^{W - 1 + e)2 



(38) 
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For our initial searches we took e = 10~^. 

Our numerical searches were carried out by minimizing a least squares likelihood function 
defined as 

6 

= ^{5vk-t\^ - Svk-kf/crl.x , (39) 
k=l 

where k indexes the six flybys reported by Anderson et al. the Svk-t^ are the theoretical 
values of the velocity discrepancies computed from our model, the Sv^-a are the observed values 
for these discrepancies reported in {2], and the ak-A are the corresponding estimated errors in these 
discrepancies given in {^]. Since the quoted ak^A values contain both systematic and statistical 
components, a least squares likelihood function is not a true statistical chi square function, but 
having a quadratic form is very convenient for the following reason. Because the theoretical values 
Svk-th are linear in the dark matter density times cross section parameters pi^e, 

6vk-th = Pi5vk-i + Pe5vk,e , (40) 

with 5vk-i^e the respective contributions from the inelastic and elastic scatterers computed with 
Pi^e = 1) the likelihood function is a positive semi-definite quadratic form in these two parameters. 
Hence for fixed values of the other six parameters ipi^e-, Ri,e, Di^e, the minimization of with 
respect to the parameters pi^e can be accomplished algebraically by solving a pair of linear equations 
in the two variables pi^e, with the result 

Pi ' 

(41) 



with coefficients given by 



Cem=2_^ 2 ' ^,m = i,e 

k=i ^fc;A 

dvk-A^Vk-/ 
- 2^ —^2 ' 

k=l ^k;A 



I, e 



(42) 

This has the effect of reducing the parameter space that must be searched numerically from an eight 
parameter space to a six parameter space, which results in a substantial saving of computational 
effort. 
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Our search procedure was then as fohows. Using a very coarse 10 point integration mesh for the 
model calculation of the flyby velocity changes, and with e = 10~^, we surveyed the six parameter 
space in 31 steps of 7r/32 for the tilt angles ipi^e, going from 7r/64 to vr — 7r/64, in 20 steps of 
2,500 km for the Gaussian centers Ri^e, going from 15,000 km to 62,500 km, and in 5 steps of 1000 
km for the Gaussian widths Di^e, going from 1,000 km to 5,000 km. For each of the 9,610,000 
steps in this survey, the values of pi^e were then optimized by using Eqs. (f4TI) and (f42]) . and the 
resulting data for values less than 25 which also had positive pe were written to a storage file. 
This left 18 potential starts for fits. For about a half dozen of these, we used the corresponding 
sets of parameter values as starting points for a six parameter minimization search using the 
CERN program Minuit, with successively 200 and then 2000 point integration meshes for the 
model calculation of the flyby velocity changes, and using double precision arithmetic throughout 
(as recommended in the Minuit documentation). Finally, using the optimized parameter values 
obtained this way, we tested for stability of the values and resulting fits with respect to program 
modifications, such as refinement of the integration mesh. The parameter space survey took several 
hours on our pentium processor laptop, the Minuit minimizations took typically minutes (or less) 
each, and the stability checks took of the order of seconds. 

This procedure showed that for e = 10^^ excellent fits could be obtained with a wide range of 
values of the radius Ri of the inelastic dark matter scatterer shell. Using the parameters for these 
good fits as a starting point, we then did a series of 5 parameter fits, each for a different fixed value 
of Ri. Also using the good fits as starting points, we did a similar series of 5 parameter fits, versus 
fixed Ri, this time with e = 10~^^ corresponding to no smoothing of the Jacobian discontinuity (up 
to the accuracy of double precision truncation errors), but using an adaptive integration program 
to adequately sample points on the trajectories where the Jacobian becomes large. These searches 
(as well as a 6 parameter fit in the e = 10~^^ case) show that the model with no smoothing has a 
distinct minimum at Ri = 34, 520 km. Results in both e cases are given in Tables II - IV. We 
caution that the e = 10~^ cases do not exactly obey the constraints between dark matter position 
and velocity required by orbital dynamics, so it is not clear at this point whether the wide range of 
Ri values and nearly exact fits obtained in this case are a reflection of just the smoothing, which 
will be present in a more realistic dark matter orbit model, or are an artifact associated with 
relaxing the orbital constraints. 

From the products piDi and PeD^ for each flt, one can use Eq. ()36p to estimate the total mass 
in the dark matter shells, in terms of the elastic and inelastic scattering parameters a^i and -Bind- 
Alternatively, given the upper bound js] on the mass of dark matter in orbit around the earth 
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between the LAGEOS satellite orbit and the moon's orbit, of 4 x IQ-^M® ~ 1.4 x lO^^GeV/c^, 
one can turn these relations into lower bounds on ad and i?inei- For example, from the values 
PiDi = 0.00304 km^ and peD^ = 19.2 km^ for fit 2d, one finds the bounds 

cTci >9.4 X lO'^^cm^ , 
Bind >l-5 X lO^^^cm^ , 

(43) 

which are consistent with the cross section range arrived at from various constraints in l|. The 

n 

spatial constraints found in pj], which require that the dark matter should be localized well away 
from the earth and the moon, are also obeyed. 

In Table V we give the results of fitting the data with Ri constrained to the value 34,520 km 
found in fit 2d (repeated in the first line of this table), versus increasing values of the Gaussian 
width Di. These results, together with those for fits 2e-g in Tables II and IV, show that the range 
of widths Di for good fits extends up to around 10,000 km. The fact that Di is not well-determined 
is also seen in the calculation leading to Table VII, where in fit 4a we give the result of repeating 
fit 2d with a refined (4000 point) integration mesh. The parameter values for fit 4a agree to within 
1% with those of fit 2d, except for Di, which in fit 4a is 2030 km, and 10^ x pi, which in fit 4a is 
1.49 km, with the product piDi matching that of fit 2d to within 1%. 

In Table VI, we show the results of basing the fit solely on a shell of inelastic scatterers, without 
a second shell of elastic scatterers. As seen, with this restriction it is not possible to get good fits, 
even when various combinations of the flyby data are excluded from the fits. For example, as shown 
on the last line of Table VI, the four parameter model with only inelastic scatterers cannot give a 
good fit to just the two flyby data points from NEAR and Messenger. In Table VII, following up 
on a suggestion by V. Toth Q], we give the results of fitting the full model, with both elastic and 
inelastic scatterers, to the flyby data, with one flyby at a time omitted from the fit. These results 
show that the predicted value for the anomaly of each omitted flyby is in qualitative accord with 
the experimental value. 

The results in Tables II - VII show that the dark matter scattering model, with inelastic and 
elastic scatterers, can account for the flyby anomaly data. One could argue that the fits are 
too good, and are indicative of "over-fitting", since there are 8 parameters in the model (9 if 
one includes e in the smoothed case), and only 6 data points. On the other hand, it was not a 
priori obvious that such a simple model should be able to account for data from a complicated 
physical process with a three-dimensional geometry, and the results shown in Table VII support 
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the view that the success of the model is not attributable to over-fitting of the data. Further steps 
in this investigation would be: (1) incorporation of further flybys into the fits, when the flyby 
parameters in Table I and the corresponding velocity discrepancy and error values are available, or 
alternatively, using fit 2d (or 4a) to predict the velocity discrepancy for future flybys, given their 
orbital parameters; (2) incorporating constraints on residual drag coming from fitting satellite drag 
measurements to conventional drag sources; (3) as suggested to me by V. Toth [gI, incorporating 
the time development of the velocity anomaly near perigree when such data becomes available from 
improved tracking of future flybys; (4) as suggested to me by J. Rosner [3], investigating possible 
constraints arising from the effect of the quadrupole moment of the dark matter shells on the 
precession of high-lying satellite orbits; (5) extending the model to include a general form of the 
weighting function it;(r, i/j); and (6) extending the model to include shells generated by precessing 
elliptical, as opposed to circular orbits, and shells generated by a precessing Schwarzschild disk 
Q]. The extensions (5) and (6), which can incorporate consistent smoothing of the Jacobian, will 
require computing resources well beyond those used here to analyze the 8 parameter model. It 
will also be necessary to address the question of mechanisms for producing dark matter shells. 
According to A. Peter the accumulation cascade suggested in pj] is not viable as a mechanism. 
Another scenario, suggested by Dr. Peter's comments and the structure of the model formulated 
here, would involve the gravitational capture by the earth of a dense (up to ~ 10^^ times galactic 
halo mean density, that is ~ 10^^ times mean ordinary matter density) condensed ball of dark 
matter into an orbit tilted with respect to earth's rotation axis; breakup of this by tidal forces 
could then lead to population of a shell of the type we have assumed.^ If the flyby anomalies are 
ultimately confirmed, detailed study of such a mechanism would be warranted. 
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APPENDIX A: CHANGES OF VARIABLE TO INTEGRATE OUT THE SPATIAL 

DELTA FUNCTION 



To eliminate the spatial delta function in Eq. (jl5p , we note that rewriting P in terms of spherical 
coordinates, 

P = r(sina;cos/3, sinwsin/?, cosw) , (Al) 
the delta function (5^(x — P) becomes 

6^{x- P) = r^^\smu;\-^S{\x\-r)S{uj{x) -io)6{(3{x) - l3) . (A2) 
Equating Eq. ()Aip with Eq. Q, we see that 



T^n :\ COS ^ COS ■ T/n ,\ ^OS 

cos^(g,-^)= sm^'(6',^) 



\/l — cos^ 9 sin^ ip Y^l — cos^ 9 sin^ ip 



(A3) 



and 

cos uj = — COS 9 sin . (A4) 

Using Eq. (jASj), on substituting Eq. ()A2p with x = x{t) into Eq. (fT5]) . we can immediately 
eliminate the r and (p integrations, leaving 

1= dt di)w{r{t),^)r{t)~'^ / d9\s\nuj\~'^ 



xF{x{t),dx{t)/dt,{GMi^/r{t))^/^U{9,4>{x{t))))6{uj{x) - uj) 



(A5) 



To carry out the 9 integration, we differentiate Eq. (|A4p . giving 

d9 —du) 



sin (jj sin ip sin 

—duj —duo 



-y/ sin^ "0(1 ~ cos^ 9) a/ sin^ V ~ cos^ u 
—rduj 



\l sin^ — z"^ 



(A6) 
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with z = rcosw. Substituting this into Eq. (|A5p . we can carry out the Q integral, leaving an the 
integral given in Eq. ()16p of the text. 

APPENDIX B: CALCULATION OF THE COEFFICIENTS C{t) AND D{t) 

From the defining cross product relations, we see that on the geocentric basis system with z 
aligned along the earth rotation axis, we have 



To express the unit velocity of Eq. (jlOp on this basis, at the intersections where x(t) = P{r, 6, cj)), 
we rewrite Eq. ([9]) as 




= {-y{t),x{t),0) 



= {- x{t)z{t),-y{t)z{t),r{tf - z{tf) 



(Bl) 



x{t)/r{t) = (cos ^ cos ^ cos — sin sine/)) 



y{t)/r{t) = (cos 9 cos ip sin (p + sin 9 cos (p) 



z{t)/r{t) 



cos 9 sin ip 



(B2) 



The third of these equations determines cos 9 and sm9 in terms of x{t) 



cos 9 



z{t) /{r{t) sin -0) 




(B3) 



while solving the first two gives sin^ and coscf) in terms of x{t) and y{t) 




(B4) 



with g = cos cos ^, h = sin^, which obey 



+ = 1- cos^ 9 sin^ 1^ = 1- z{tf/r{tf 



(B5) 
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Substituting Eqs. ()B3P and ()B4p into ([TO]) gives the velocity components at the intersections 
expressed in terms of x{t), and comparing with Eq. (jBlh then identifies the coefficients C{t) and 
D{t) appearing in the decomposition of U± on the intrinsicahy defined basis ny and h±. 
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TABLE II: Flyby anomaly fits 







GLL-I 


GLL-II 


NEAR 


Cassini 


Rosctta 


Messenger 


5vA (mm/s) 




3.92 


-4.6 


13.46 


-2 


1.80 


0.02 


a A (mm/s) 




0.3 


1.0 


0.01 


1 


0.03 


0.01 


6vth fits la-e 


< 10-6 


3.92 


-4.60 


13.46 


-2.00 


1.80 


0.020 


(5wth fit 2a 


2.07 


3.98 


-5.5 


13.46 


-3.1 


1.79 


0.021 


Svth fit 2b 


1.68 


4.15 


-5.2 


13.46 


-2.9 


1.80 


0.020 


Svth fit 2c 


1.29 


4.13 


-5.0 


13.46 


-2.8 


1.80 


0.020 


(5wth fit 2d 


0.51 


3.90 


-4.6 


13.46 


-2.7 


1.80 


0.020 


Svth fit 2e 


0.52 


3.88 


-4.6 


13.46 


-2.7 


1.80 


0.020 


(Jvth fit 2f 


0.70 


3.84 


-4.7 


13.46 


-2.7 


1.80 


0.021 


5vth fit 2g 


7.5 


3.76 


-4.7 


13.46 


-2.8 


1.73 


0.028 



Fits la-e are for the smoothed model with e = 10"^ and trapezoidal integration, resulting from a 
five-parameter fit with Ri constrained to the values shown in Table III. Fits 2a-g are for the un-smoothed 
model (e = 10~^^, which is below truncation errors) and adaptive trapezoidal integration, resulting from a 
five-parameter fit with Ri constrained to the values shown in Table IV. 



TABLE III: Parameter values for fits la-e 



fit 


10*3 X Pi (km) 


102 


X pe (km) 


ipi (rad) 


tjje (rad) 


Ri (km) 


Di (km) 


Re (km) 


De (km) 


la 


0.304 




0.268 


1.926 


0.3939 


30000 


6278 


28620 


6303 


lb 


1.55 




0.245 


1.261 


0.3945 


40000 


2185 


27985 


5890 


Ic 


0.411 




0.261 


1.374 


0.3952 


50000 


13540 


28450 


6299 


Id 


0.351 




0.253 


1.381 


0.3946 


60000 


20193 


28340 


6334 


le 


0.343 




0.248 


1.394 


0.3942 


70000 


25780 


28240 


6367 



TABLE IV: Parameter values for fits 2a-g 



fit 


10^ X Pi (km) 


10^ X Pe (km) 


i^i (rad) 


V'e (rad) 


Ri (km) 


A (km) 


Re (km) 


De (km) 


2a 


0.537 


0.323 


1.767 


0.3902 


25000 


3030 


29370 


6678 


2b 


0.827 


0.316 


1.626 


0.3902 


30000 


3030 


29370 


6678 


2c 


0.965 


0.309 


1.515 


0.3902 


32500 


3030 


29370 


6678 


2d 


1.000 


0.288 


1.372 


0.3902 


34520 


3030 


29370 


6678 


2c 


0.655 


0.288 


1.369 


0.3902 


35000 


4663 


29370 


6678 


2f 


0.348 


0.288 


1.364 


0.3902 


37500 


9223 


29370 


6678 


2g 


0.290 


0.286 


1.361 


0.3902 


40000 


11681 


29370 


6678 
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TABLE V: Flyby anomaly fits with Ri = 34520 and indicated values of Di 







GLL-I 


GLL-II 


NEAR 


Cassini 


Rosctta 


Messenger 


Sva (mm/s) 




3.92 


-4.6 


13.46 


-2 


1.80 


0.02 


dvth Di = 3030 


0.51 


3.90 


-4.8 


13.46 


-2.7 


1.80 


0.02 


Svth D, = 6060 


0.68 


3.94 


-4.8 


13.46 


-2.8 


1.80 


0.02 


(5wth Di = 9090 


1.3 


3.97 


-5.1 


13.46 


-3.0 


1.80 


0.02 


(5i;th Di = 12120 


4.2 


3.88 


-4.7 


13.46 


-4.0 


1.80 


0.02 



TABLE VI: Flyby anomaly attempted fits with only inelastic scatterers 







GLL-I 


GLL-II 


NEAR 


Cassini 


Rosctta 


Messenger 


5vA (mm/s) 




3.92 


-4.6 


13.46 


-2 


1.80 


0.02 


(5wth fit 3a 


0.63 X 10-^ 


1.87 


1.8 


13.0 


1.5 


2.9 


2.4 


(5wth fit 3b 


0.63 X 10'"^ 


1.87 




13.0 




2.9 


2.4 


5vt\i fit 3c 


0.16 X lO'' 


1.93 




13.4 




3.0 




5vt\i fit 3d 


0.61 X 10^ 






13.0 






2.5 



Fit attempts with only inelastic scattering; entries labeled - were excluded from the corresponding fit. 
To two decimal places, all fits in this table correspond to the parameter values p, = 0.14, ■^j = 1.13, 
Ri = 40000, and A = 2000. 



TABLE VII: Flyby anomaly fits to five of the six flybys 





\' 


GLL-I 


GLL-II 


NEAR 


( 'asriiiii 




M(>riS('llg(T 


5va (mm/'s) 




3.92 


-4.0 


i:-!.4() 


_2 


1.80 


0.02 


^Wth fit 4a 


0.49 


3.90 


-4.6 


13.46 


-2.7 


1.80 


0.02 


5vth fit 4b 


0.45 


3.71 


-4.4 


13.46 


-2.6 


1.80 


0.02 


5v\,h fit 4c 


0.49 


3.91 


-4.6 


13.46 


-2.7 


1.80 


0.02 


5vth. fit 4d 


0.40 


3.93 


-4.4 


16.03 


-2.6 


1.80 


0.02 


5vih. fit 4e 


0.63 X 10-3 


3.92 


-4.6 


13.46 


-2.7 


1.80 


0.02 


(Jwth fit 4f 


0.40 X 10-1 


3.93 


-4.5 


13.46 


-2.2 


1.62 


0.02 


5vi\, lit 4g 


0.19 


;-').94 


4.:-! 


i:-!.4() 


-2.3 


1.80 


0.12 



Fit 4a is a fit with all six flybys included. Fits 4b-4g are fits with one flyby at a time excluded; the 
predicted value for the flyby omitted in each fit is in boldface. These fits use a factor of 2 finer mesh than 
fit 2d. 



